function [J,g_norm,w,delta_w,Solution,u_max]=newton(U,delta_U,Solution)
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

[g_t,delta_g_t,g_p,g_norm,res,pu,delta_pu]=gradient(1,U,delta_U,1);
[N_D,delta_N_D,n_d]=newton_direction(U,delta_U);
w=U;
delta_w=delta_U;
g_norm_s=g_norm;
epsilon=input('Starting Newton Method. Input its accuracy = ');
it_max=input('Input Maximum Iteration Number it_max=');
disp('     it_n      g_norm    err_max   energy    u_max at       (x,y)           s');
it_n=1;
while g_norm>epsilon
    s=2.0;
    while g_norm_s-g_norm >= -0.5*s*g_norm;
        w_s=w-s*N_D';
        delta_w_s=delta_w-s*delta_N_D';
        [g_s_t,delta_g_s_t,g_s_p,g_norm_s,res,pu,delta_pu]=gradient(1,w_s,delta_w_s,1);
        s=s/2;
    end
    J=-energy(1,w_s,delta_w_s,1);
    err_max=max(abs(res));
    [u_max,II]=max(abs(w_s));
    disp([it_n g_norm_s err_max J u_max x_t(II) y_t(II) s]);
    it_n=it_n+1;
    w=w_s;
    delta_w=delta_w_s;
    Solution=Solution-2*s*n_d;
    g_norm=g_norm_s;

    u_max=max(Solution);
    u_min=min(Solution);
    set(gca,'ZLim',[u_min u_max]);
    %pdesurf(p,t,Solution); drawnow;
    pdeplot(p,[],t,'XYData',Solution,'ZData',-Solution,'ColorBar','off'); drawnow;
    if it_n-1>it_max
        disp('It_n > It_max');
        break;
    end
    [N_D,delta_N_D,n_d]=newton_direction(w,delta_w);
end
J=energy(1,U,delta_U,1);
J=-J;
